# 
# Description:
# This script constructs Table A3
# for Appendix A
#
# Requires: pres_all_terror.dta, CAINC30__ALL_STATES_1969_2017.csv, county_ucr_offenses_known_yearly_1960_2017.rda, terror.dta
#
# Author: Sean Nossek
#
# Date: August 04, 2020

# install pacman
if (!require("pacman")) install.packages("pacman")

#lfe version manage
packageurl <- "https://cran.r-project.org/src/contrib/Archive/lfe/lfe_2.8.tar.gz"
install.packages(packageurl, repos=NULL, type="source")
library(lfe)


# load other packages
pacman::p_load(tidyverse,
               readstata13,
               noncensus,
               xtable,
               stargazer)

# set wd
setwd("C:/Users/Sean/Dropbox/Eran/Terrorism_Left/Paper/Revisions/Revised Paper/Replication/Paper Data")

# create state_capital dummy
pres_terror <- read.dta13("pres_all_terror.dta") %>%
  mutate(fips = as.numeric(paste0(state_fips, str_pad(county_fips, 3, pad = "0"))))

# region data
data(states)

states <- states %>%
  dplyr::select(state, region) %>%
  rename(state_abbrev = state)

# merge regions
pres_terror <- merge(pres_terror, states, by = "state_abbrev", all.x = TRUE) %>%
  mutate(GeoFIPS = str_pad(fips, 5, pad = "0"))

#load bea county level data
bea <- read.csv("CAINC30__ALL_STATES_1969_2017.csv", stringsAsFactors = FALSE) %>%
  gather(year, val, X1969:X2017) %>% 
  mutate(Description = paste(Description, Unit, sep = " - ")) %>%
  dplyr::select(-c(LineCode, Unit)) %>% spread(Description, val) %>%
  dplyr::select(-c(` - `, TableName, IndustryClassification)) %>%
  mutate(GeoFIPS = ifelse(str_length(GeoFIPS) != 6, NA, str_remove(GeoFIPS, " ")), year = as.numeric(str_remove(year, "X"))) %>%
  ungroup()

#select region
static <- pres_terror %>% 
  ungroup() %>%
  dplyr::select(fips, GeoFIPS, region) %>%
  distinct()

#merge region with bea data to have county-year panel
balance <- merge(bea, static, by = c("GeoFIPS")) %>%
  arrange(fips, year)

# read crime data
load("county_ucr_offenses_known_yearly_1960_2017.rda")

crime <- county_ucr_offenses_known_yearly_1960_2017 %>% 
  ungroup() %>%
  rename(GeoFIPS = fips_state_county)

# merge this
balance <- merge(balance, crime, by=c("GeoFIPS","year")) %>% 
  filter(year <= 2013)

# add which counties had attacks, how many were successful, for balance tests
success_counties <- read.dta13("terror.dta") %>%
  group_by(FIPSCounty, year) %>%
  summarise(success = sum(success)) %>%
  mutate(GeoFIPS = str_pad(FIPSCounty, 5, side = "left", "0"), attack = 1)

balance <- merge(balance, success_counties, by = c("GeoFIPS","year"), all.x = TRUE) %>%
  mutate(success = ifelse(is.na(success), 0, success),
         attack = ifelse(is.na(attack), 0, attack)) %>% 
  group_by(GeoFIPS) %>%
  mutate(success_lead = lead(success, order_by = year),
         column_a = as.numeric(success_lead >= 1), 
         column_b = as.numeric(attack == 0 & column_a == 0),
         success_test = ifelse(column_a == 1, 1, ifelse(column_b == 1, 0, NA)),
         keep_1969 = ifelse(year != 1969, 1, ifelse(success_test == 1, 1, 0))) %>%
  filter(keep_1969 != 0)

# clean variables for comparison
balance <- balance %>% mutate(l_jobs_per_cap = log(as.numeric(`Total employment (number of jobs) - Number of jobs`) / as.numeric(` Population (persons) 3/ - Number of persons`)),
                              l_earnings_per_cap = log((as.numeric(`  Wages and salaries - Thousands of dollars`)*1000) / as.numeric(` Population (persons) 3/ - Number of persons`)),
                              l_retirement_per_cap = log(as.numeric(`  Per capita retirement and other 4/ - Dollars`)), l_unemp_insurance_per_cap = log(as.numeric(`  Per capita unemployment insurance compensation 4/ - Dollars`) + .001),
                              l_population = log(as.numeric(` Population (persons) 3/ - Number of persons`)), l_officers_killed_felony_per_cap = log( officers_killed_by_felony  + .001/ as.numeric(` Population (persons) 3/ - Number of persons`)), l_murder_per_cap = log( actual_murder +.001/ as.numeric(` Population (persons) 3/ - Number of persons`)),
                              l_mtr_veh_theft_per_cap = log( actual_mtr_veh_theft_total +.001/ as.numeric(` Population (persons) 3/ - Number of persons`)), l_robbery_per_cap = log( actual_robbery_total +.001/ as.numeric(` Population (persons) 3/ - Number of persons`)),
                              l_violent_crime_per_cap = log( actual_index_violent +.001/ as.numeric(` Population (persons) 3/ - Number of persons`)), l_prop_crime_per_cap = log( actual_index_property +.001/ as.numeric(` Population (persons) 3/ - Number of persons`)),
                              l_crime_per_cap = log( actual_index_total +.001/ as.numeric(` Population (persons) 3/ - Number of persons`)), northeast = as.numeric(region == "Northeast"), south = as.numeric(region == "South"), midwest = as.numeric(region == "Midwest"), west = as.numeric(region == "West"))

# save controls for use in robustness checks

saveRDS(balance, "county_balance.RDS")

# make success v failure dummy
balance <- balance %>%
  mutate(success_v_fail = ifelse(success >= 1, 1, ifelse(attack == 1, 0, NA)))


##############################
######### OLS FE Regression ##
##############################

varlist <- c("l_jobs_per_cap", "l_earnings_per_cap", "l_population", "l_retirement_per_cap","l_unemp_insurance_per_cap", "l_officers_killed_felony_per_cap", "l_murder_per_cap","l_mtr_veh_theft_per_cap","l_robbery_per_cap","l_violent_crime_per_cap","l_prop_crime_per_cap","l_crime_per_cap")

fail_test <- lapply(varlist, function(x) {
  felm(as.formula(substitute(success_v_fail ~ i | GeoFIPS + year | 0 | GeoFIPS, list(i = as.name(x)))), data = balance)
})


# Table A3
stargazer(fail_test)

